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COMPUTATION OF SYNTHETIC SPECTRA FROM SIMULATIONS OF RELATIVISTIC SHOCKS 
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ABSTRACT 

Particle-in-cell (PIC) simulations of relativistic shocks are in principle capable of predicting the 
spectra of photons that are radiated incoherently by the accelerated particles. The most direct method 
evaluates the spectrum using the fields given by the Licnard-Wicchart potentials. However, for 
relativistic particles this procedure is computationally expensive. Here we present an alternative 
method, that uses the concept of the photon formation length. The algorithm is suitable for evaluating 
spectra both from particles moving in a specific realization of a turbulent electromagnetic field, or 
from trajectories given as a finite, discrete time series by a PIC simulation. The main advantage of 
the method is that it identifies the intrinsic spectral features, and filters out those that are artifacts 
of the limited time resolution and finite duration of input trajectories. 

Subject headings: radiation mechanisms: general — methods: numerical — gamma-ray burst: general 



relativistic processes 



1. INTRODUCTION 

Particle acceleration at relativistic shocks is thought 
to be responsible for the high-energy nonthermal pho- 
tons observed from a variety of astrophysical objects, 
such as gamma-ray bursts, pulsars and blazars. To test 
this hypothesis, reliable predictions of the photon spec- 
tra are needed. Analytic models have been used to pro- 
vide estimates of the expected asymptotic power-law in- 
dex at high energy and the maximum attaina ble photon 
energy (|Derishevi 120071 iKirk k Reville! 1201 Of ), but they 
cannot currently take account of potentially important 
effects, such as the role of self-generated turbulence in 
the vicinity of the shock. Particle in cell (PIC) simula- 
tions, on the other hand, have the potential to capture 
these effects, and have recently begun to provide evi- 
dence that Fermi acceleration is a natural consequen ce 
of relativistic shock formation dSpitkovskv 2005, 2 008al lbl; 
iSironi k Spitkovskvl |2009al IMartins et al.l l2009bf) . In 
principle, these simulations are capable of reproducing 
the essential physics; they are ab initio in the sense that 
all processes are reproduced by evolving the electromag- 
netic fields and the particle distribution according to the 
classical equations of motion and Maxwell's equations. 

However, to compare the simulations with observa- 
tions, it is essential to understand the predicted ra- 
diative signatures. Using results from PIC simula- 
tions, several grou ps have comp uted the emission at 
relativistic shocks (IHededall 120051 IMartins et al.ll2009at 
ISironi k Spitkovskvl l2009bHMedvedev et al.l 120101) . but 
the results show substantial differences. This is not nec- 
essarily due to the way in which the spectra were evalu- 
ated, since the electromagnetic fields and energetic par- 
ticle distribution vary strongly from simulation to sim- 
ulation. On the other hand, it does not rule out such 
a dependence. In each case, the method employed to 
compute the emission is the same: the electric field pro- 
duced at a virtual detector by a single particle trajectory 
is evaluated using the Lienard-Wiechart potentials, the 
result is Fourier transformed and then averaged over a 

brian.reville@mpi-hd.mpg.de, john.kirk@mpi-hd.mpg.de 



large family of trajectories. This procedure is expensive 
in terms of computing resources, especially if one wants 
to compute the high-energy emission of relativistic par- 
ticles, because of (i) the extremely high time resolution 
required to describe high-energy photons, (ii) the large 
number of virtual detectors required to resolve the nar- 
row radiation beam of a relativistic particle (which scales 
as 7 2 ) and (iii) the long time series needed to account for 
low-frequency emission. 

In this paper, we present an alternative approach. 
For a given observing frequency, we identify at each 
point on a particle trajectory the length that con- 
tributes coherently to the emission. In a quantum pic- 
ture, this is known as t he photon formation length 
(jAkh lezer k ShurgallT987l) . Using this as a guide, we 
then perform the integrations using a new algorithm that 
is optimized for highly relativistic particles. 

The essential information on the particle trajectory 
can be supplied to the algorithm in two different ways: 
If the electromagnetic fields are prescribed as a func- 
tion of space and time, then the trajectory can be in- 
tegrated using standard adjustable-step methods. As 
examples, we present in section [3] computations of the 
emission spectrum from isotropic particle distributions 
immersed in stationary, turbulent magnetic fields. In 
this case, the algorithm is employed in the inner loop of 
a multi-dimensional integration, which is performed us- 
ing a Monte-Carlo method. On the other hand, if the 
fields are not known at all points in space and time, in- 
terpolation is required. This is the case, for example, 
in PIC simulations, where the fields, particle positions 
and velocities are known only at discrete times and loca- 
tions. We discuss this situation and suggest a procedure 
for implementing the algorithm in section |4j 



2. EQUATIONS FOR THE EMISSIVITY 

In the classical theory of electrodynamics, the spectral 
and angular distribution of radiation produced by a sin- 
gle particle in vacuum in the direction n is given by the 
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well-known formula (e.g. lLandau fc LifshitdU971[ ) 



dE 



dcudtt 4ir 2 c 



<r 



+00 n x 



(n - 0) x 



^ifcjt — k-x 



(l-n-/3) 2 



W)d<l) 



where k = um/c. Equation ([1} is usually used as the 
starting point for the numerical computation of radiation 
signatures from PIC cod es — a de t ailed description of the 
method can be found in lHeded"al (|2005[ ). However, there 
are three disadvantages of this form of the emissivity 

1. the term e 1 "* is rapidly oscillating 

2. the term 1 — n ■ f3(t) in the denominator produces 
a very sharply peaked function when used for the 
trajectory of a relativistic particle 

3. the range of integration extends over the entire sec- 
tion of the trajectory on which the acceleration is 
nonzero, making it difficult to relate the expres- 
sion to a local emissivity and, hence, to compute 
time-dependent emission. 

Straightforward transformations lead to a number of 
alternative forms for Equation ([lj, for example, 
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P(n,w,t)dt 



dr[l-f3(t + r)-/3(t)} 



(2) 



cos(u}[r-n-(x{t + r)-x(t))/c]) .(3) 

This expression is exact, and has the advantage that, 
provided variations on the timescale a; -1 are small, the 
quantity P(n,ui,t), when suitably averaged, can be in- 
terpreted as the instantaneous spectral pow er radiated 
per u nit solid angle about the direction n (iSchwingerl 

MM- 

In practice, it is necessary to truncate the integrals in 
and (13]) to finite intervals. From the form of the in- 
tegrand in (|3|), it is clear that the endpoints should be 
chosen such that at least the first few periods of the co- 
sine function are included. This leads to the concept of 
the photon formation time or coherence time, which ap- 
plies to both the quantum a nd classical formulations o f 
the problem (for a review, see lAkhiezer fc Shurgalll987l ). 
For a given Fourier mode, with wavelength A = 2ttc/lu, a 
particle trajectory contributes coherently to the instan- 
taneous power radiated at time t until it has lagged at 
least one wavelength behind the wavefront emitted at 
time t. Thus, the coherence or formation time r co h is 
determined implicitly by the equation 



w (r coh - \x(t + T coh ) - x(t)\ /c) = 2tt 



(4) 



To compute the radiated power, one needs to know the 
trajectory accurately over several coherence times. For 
relativistic particles, a wavefront can take a consider- 
able amount of time to separate one wavelength from 
the particle, particularly at low frequencies, when the 
wavelength is long. 



In the context of Fermi acceleration at relativistic 
shocks, an angular dependent calculation of the emis- 
sion from an individual particle is unnecessary, provided 
one is interested only in the high-energy emission from 
accelerated particles. This is because the characteris- 
tic radiation beaming angle of I/7 is much smaller for 
these particles than the scales on which anisotropy in the 
particle distribution can be expected, which is roughly 
the reciprocal o f the Lorentz factor of the fluid motion 
into the shock (jAchterberg et al] 120011 ). Hence, when 
summed over all plasma particles, these sharp emission 
peaks are smoothed out. In this case, it is advantageous 
to work with an angle-integrated expression for the indi- 
vidual particle spectrum. Integrating equation ([3]) over 
solid angle, gives: 
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where 
P(u,t) 
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with 



F(.,t,r) = Sin[ " (T - A) ]- Sin[ " (r + A)] (7) 



A(t,r) = \x(t + T)-x(t)\/c 



(8) 



Equation ^ is also exact, and P(uj, t) can be interpreted 
as the power radiated at time t in unit angular frequency 
range, again subject to the condition th at it varies slowly 
on the timescale w _1 (jSchwingerl 119491 ) . This condition 
is not always fulfilled for the trajectories we consider. 
In particular, it is violated when the acceleration felt by 
the particle fluctuates rapidly whilst the velocity remains 
within the beaming angle of the radiation ("jitter" radi- 
ation). Nevertheless, equation §5§ for the total radiated 
energy remains valid, although P(oj,t), which we call the 
"instantaneous power" , is not necessarily positive defi- 
nite. 

2.1. Computation of the instantaneous power 

For relativistic particles, and for small r, such that the 
particle displacement A defined in (JSj) is approximately 
/3|r|, the function F(uj,t,r) in (0 contains two kinds 
of term: those that oscillate rapidly in r with frequency 
~ to, and those that oscillate slowly, with frequency w/7 2 . 
Physically, the latter arise because the particle chases the 
wavefront, remaining close to it for a relatively long time. 
It is convenient to separate these terms: 

P{u,t) = P 1 (u,t) + P 2 {u,t) (9) 
Pi (w, t) = ^ J°° dr [1 - f3(t) ■ f3(t + t)} 
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In the case of P2, the time r over which the trajectory 
is sampled is very short, ~ 1/w, whereas in the case of 
Pi it is much longer, ~ 7 2 / w - To develop an approxi- 
mation scheme for these terms we introduce quantities 
that describe the deviation of the trajectory from ballis- 
tic motion. These are the deviation in position 

8x(t,T) = x(t + T)-x(t)-CT/3(t), (12) 

the deviation in velocity: 

5f3(t, T )=(3(t + T)-f3(t) (13) 

and the deviation of the displacement 

5A(t,T)=A(t,T)-\r\p(t)- (14) 

Clearly, to zeroth order in these deviations, the instan- 
taneous power must vanish, since a particle undergoing 
uniform motion does not radiate. Furthermore, because 
of the relatively long sampling time, the dominant higher 
order contributions in the deviations come from Pi . For 
frequencies large compared to the instantaneous angu- 
lar frequency (the local gyrofrequency) , the higher order 
contributions in the P2 term can be neglected to give 
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where we have introduced the phase-lag g(u>,t,T): 
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At this point it would be possible to proceed by evalu- 
ating analytically the integral involving the second term 
in ([TBI: 
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Indeed, iSchwingen () 19491 ) followed this path in deriv- 
ing an analytic expression for the synchrotron emissivity. 
However, because the integrands are oscillatory, it is in- 
stead preferable to group them together. Transforming 
the integration variable from r to the phase-lag g de- 
fined in (|16p . in the case of the first term in (|15| . and as 
g = ljt(1 + (3(t)) in the case of the second term, leads to 
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As required, P vanishes to zeroth order in the deviations 
from a ballistic orbit, (O, (|T3J) and ([Tl]). The grouping 
of the terms in Equation f| 18|) in this manner is especially 
important at high frequencies, where the higher order 
terms in Pi and P2 are small. In this limit, the two 
terms can be expressed as 



lim Pi, 2 = ± : 



2c i{tym 



(21) 



and cancel exactly when summed. In a numerical eval- 
uation, a small error remains, which grows linearly with 
bj. Grouping the terms together prevents the growth of 
this error. 

Under the assumptions that the electromagnetic fields 
vary slowly on the timescale of a photon forma- 
tion length, and that linear acceleration emission (e.g. 
lSchwingerlll949[) is unimportant, we demonstrate in ap- 
pendix [A] that (fl8|) reduces to a local emissivity. This is 
an obvious generalization of standard synchrotron emis- 
sion, which takes account of acceleration in both mag- 
netic and electric fields by formulating it in terms of the 
local curvature of the trajectory: 



P(uJ,t)- 



2tt uj c 



where 



u) c = 3j 3 ck/2 



dxK 5/3 (x) (22) 



(23) 



and the curvature k is defined locally in terms of the 
particle velocity and acceleration j3 and 0: 



0x0 



cf3 3 



(24) 



A perturbative approach that includes linear accelera- 
tion emission as a firs t orde r correction to (|2"2"|) has been 
presented bv lMelrosel (| 19781 ) . 

To perform the integration in f| 18|) numerically, we first 
split it at the points where sing = 0, i.e., g = nn, (n = 
0,±1,±2...), and write it as an infinite sum 
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and t and g arc considered to be functions of g and t, 
defined implicitly in (fT"6|) and (|20|) . According to its defi- 
nition (U, integration from g — to g — 2-7T corresponds 
precisely to integration over one photon formation time. 
We therefore anticipate on physical grounds that taking 
the first few terms should give a good approximation. 
However, the function Q(g,t) grows linearly with g for 
small g, before decreasing monotonically above some crit- 
ical value g* . In this case, the sum in (|25|) does not begin 
to converge until n > n* = g* /n. For constant cur- 
vature, it is straightforward to show that g* w 3w/u; c , 
so that n* becomes large only if one tries to compute 
the emissivity well above the cut-off frequency. In gen- 
eral, we have found that a substantial improvement can 
be achieved b y employing the E uler-van Wijngaarden 
transform (e.g. iPress et al.lll986l section 5.1) to acceler- 
ate the convergence, whilst retaining a minimum number 
of about 20 terms in order to preserve accuracy at high 
frequencies, where the power radiated is low. 

Evaluation of the instantaneous power based on (|25j) 
requires knowledge of the functions j3(r) and 6A(t). In 
the next section we apply this approach to finding the an- 
gular integrated emission of an isotropic, mono-energetic 
particle distribution in prescribed, stationary, turbulent 
fields, in which these functions can be found using an 
adjustable-step integration of the trajectory. In section^ 
on the other hand, we discuss the application to a tra- 
jectory that is known only as a discrete time series, for 
example, a trajectory from a PIC simulation. 

3. PRESCRIBED FIELDS 

3.1. Isotropic particle distribution 

Equation ([5|) describes the energy emitted by a single 
particle. If we now consider the possibility of ./V parti- 
cles emitting incoherently whilst following trajectories in 
a prescribed field in a volume V, and allow them to do 
so for a time T, then the average power L emitted by 
these particles is obtained by summing over the individ- 
ual contributions: 
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where Pi(t) is the instantaneous power of the i'th par- 
ticle. Replacing the sum by an integral over the ex- 
act (Klimontovich) phase space distribution /k(x,p, t) = 

Z)j=i 5[x- Xi(t)] S\p- Pi(t)] where Xi(t),pi(t) are the 
phase-space coordinates of the i'th particle at time t, 
leads to 
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where P{xi(t), Pi(t), t) = P t (t). 

In general, both the electromagnetic fields that deter- 
mine the particle trajectories and the phase space distri- 
bution fluctuate in time. However, L is a time-averaged 
quantity. If we are interested in the emission from a sys- 
tem containing prescribed, static fields, then P(x, p, t) is 
not an explicit function of time, so that 
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j A 3 xA 3 pP(x,p)± j T ^Atf K (x,p,t) (29) 



If, in addition, we look at the radiation from a station- 
ary coarse-grained particle distribution f(x,p), then, re- 
placing the time-averaged Klimontovich function by this 
distribution leads to 
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(30) 



For the case of fluctuations in only the magnetic field, for 
example, the particle energy is an integral of motion, and 
any homogeneous, isotropic function of the Lorentz fac- 
tor 7(p) is a stationary solution of the kinetic equation. 

Setting f(x,p) = $^£(7 ~ 7(p)), wc nn d 
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where fl = p/p. Thus, in order to compute the power 
radiated per unit frequency interval, we must integrate 
the instantaneous power over all directions of the veloc- 
ity vector at each point and over all positions within 
the source. In the following subsections we present com- 
putations of the radiation produced from an ensemble 
of relativistic particles in static turbulent magnetic field 
configurations, employing a Monte Carlo integration of 
Equation ([3i~jl . 



3.2. Emission spectrum 

The character of the radiation produced by a relativis- 
tic particle depends on whether the strength parameter 



eFX 



(32) 



is greater than or less than unity, where A is the typical 
size of the field structures and eF = (dp \ /At) the aver- 
age t ransverse force on the particle (jLandau fc Lifshita 
Il971h . This Lorentz invariant parameter is analogous to 
the strength parameter commonly used in laser plasma 
physics, and is also sometimes called the "wiggler" or 
"undulator" parameter. For static fields, it determines 
roughly the ratio of the deflection angle to the beaming 
angle for a particle traversing a typical structure. For 
simplicity, electric fields are neglected for the remainder 
of this section (F = B±). Typically, the magnitude of 
the strength parameter determines whether the particle 
radiates in the synchrotron regime (a > 1) or in the so- 
called jitter/diffuse synchrotron regime (a < 1). 

For a given B± and A, the maximum photon energy can 
be determined. However, the full details of the spectrum 
produced by a particle, even in a relatively simple field 
configuration, can be quite complicated. Using the algo- 
rithm presented in the previous section, the equations of 
motion can be integrated simultaneously with equation 
(|25p . providing the complete spectrum. For the results 
that follo w a fifth-order ada ptive Runge-Kutta integrator 
was used ([Press et al.lll986[ ). With the aid of some illus- 
trative examples, we demonstrate how different spectral 
features can be produced, and emphasize the properties 
of the fields required to do so. 

3.3. Uniform fields - the synchrotron approximation 

As a first example, the radiation produced from a parti- 
cle gyrating in a uniform field is compared to the analytic 
solution for synchrotron radiation, Eq. (|22| . The results 
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Fig. 1. — Instantaneous power spectrum produced by a particle 
of Lorentz factor 7 = 10 2 with different pitch angles in a uniform 
magnetic field B, as a function of angular frequency in units of 
uj g = eB/mc, with 9 the angle between the particle velocity and the 
field. The numerical values (solid lines) are in excellent agreement 
with the instantaneous synchrotron approximation (dashed lines). 
The noise at high frequencies can be controlled by increasing both 
the number of terms taken in the series and the integrator accuracy. 

are shown in Fig. [T] and are in excellent agreement with 
the analytic result. At very high frequencies, the forma- 
tion lengths become extremely short, and Q can be linear 
in g for several periods. Errors in the particle integrator 
can also become an issue. For the results shown in Fig. [T] 
a fractional error control of 10 -7 was used and the power 
was summed from n = 0ton = ±10. As discussed in sec- 
tion 12.1 [ the series fails to converge if n* > 10. However, 
this occurs only well above the cut-off frequency. As we 
discuss in section |4j when dealing with discrete time se- 
ries, the "synchrotron approximation" must be taken at 
high frequencies, where the formation length is small. 

The grouping of the terms described in Equation (|18p 
is vital in keeping the high frequency noise below the in- 
tegrator accuracy. This is also important for calculations 
in turbulent fields when there is a large variation in the 
roll-over frequency of the instantaneous power. Provided 
the high frequency noise remains below the threshold, the 
results are reliable. 

3.4. Turbulent fields 

In a turbulent magnetic field, the particle trajectory 
and r esulting radiation spectrum ar e generally quite com- 
plex (|Toptvgin fc Fleishman! H987I ). Nevertheless, sev- 
eral qualitative features can be understood in terms of 
the strength parameter, although the product BX is re- 
placed by a different value for each Fourier mode. There 
is now no single strength parameter but rather a spec- 
trum a(k) = 2ireB(k)/ (mc 2 fc), and radiation produced 
depends on several factors, most notably the turbulent 
spectrum and the magnitudes of a(fc m i n ) and a(fc max ). 

To investigate the effect of different turbulent field 
parameters, static fields are constructed with the re- 
quired properties. This is done using a discrete 
Fourier transform descrip tion following the method of 
iGiacalone fc Jokipiil (I1999D . The magnetic field at a po- 
sition x is B(x) = Bo + SB(x), where -Bo represents an 
external uniform mean field. The turbulent field com- 
ponent is generated using N Fourier modes, each with a 



random phase, direction and polarization. In the limit 
of large N, 



N 



SB(x) = lim 



71=1 



(33) 



represents an isotropic turbulent field. Here A n , j3 n , k n 
and £ n are the amplitude, phase, wave vector and polar- 
ization vector for each mode n respectively. The polariza- 
tion vector is determined by a single angle < ip n < 2n 



cosip n e x + ismip n e y 



(34) 



where e x and e y are vectors, orthonormal to e z = k n /k n . 
The vector k n is determined by two additional angles, 
< 8 n < 7r and < <p n < 2n, and, for an isotropic 
distribution, should be uniformly distributed on the unit 
sphere. These two angles define a rotation matrix t hat 
determines e x and e y (e.g. IGiacalone fc Jokipiilll999|) . 
The amplitude of each mode is 



A 2 n = a 2 G n 



N 



where the variance a 2 is chosen such that the turbulent 
field is normalized to give the required turbulence level: 



(SB 2 ) 



B 2 



(SB 2 ) 



(35) 



We use the following form for the power spectrum 



G n — 



1 + (k n L c ) a 



(36) 



where L c is the correlation length of the field and a is the 
asymptotic spectral index of the turbulence spectrum. 
For the three-dimensional fields used in this paper the 
normalization factor is AV n = 47r/c 2 Afc„, and the Ak n 
are chosen such that there is equal spacing in logarith- 
mic fc-space, over the finite interval fc m in < k < /c max - 
For a detailed discussion of the stat istical properties o f 
fields constructed in this manner see ICasse et all (120(51) . 
The field can be constructed at any point in space by 
summing over the N modes, providing an infinite spa- 
tial description of the fields. This avoids the need for 
boundary conditions. The parameters used for each field 
construction are given in Table [TJ 

The spectra are produced using a Monte Carlo inte- 
gration of Equation (|31[) . At each frequency u> a sample 
particle of fixed Lorentz factor is placed at a random lo- 
cation Xi inside a volume with dimensions several times 
the size of the correlation length, L c , of the turbulent 
field. To represent an isotropic particle distribution, the 
particle is given a random initial direction f2j, and the 
instantaneous power Pj is calculated. The average power 
emitted per particle at each frequency is determined us- 
ing a Monte Carlo integration: 
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TABLE 1 
Turbulent Field Parameters 



Field 


-^rms 


2tt/ &max 


27r//c m ; n 


L c 


V 


a 


A 


1.0 


2 


160 


80 


1.0 


11/3 


B 


1.0 


2 


320 


160 


1.0 


11/3 


C 


0.05 


0.5 


10 


5 


1.0 


11/3 


D 


1.0 


0.05 


10 


5 


1.0 


8/3 


E 


0.1 


0.1 


1 


0.5 


1.0 


9/3 


F 


1.0 


0.05 


10 


5 


0.1 


8/3 


G 


1.0 


0.1 


10 


5 


0.9 


8/3 



Note. — Parameters used in the field constructions for 
turbulent field spectra. All quantities are dimensionless, with 
the magnetic field in units of an arbitrary normalization value 
Bq. All length scales are in units of mc 2 /eBq. The maximum 
strength parameter in each run is given approximately by the 
product 27rB rms /fc min . 
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Fig. 2. — Radiation spectra for an isotropic distribution of mono- 
energetic particles in a fully turbulent field with strength parameter 
a S> 1, as a function of frequency in units of u) s = eBo/mc. dW/dui 
represents the average power emitted by each particle and is plot- 
ted in units e 2 u) g /2nc. The dashed lines show the integral of the 
instantaneous power, evaluated in the synchrotron approximation. 

The number of integration points n is determined from 
the condition that the standard deviation error estimate 
is well below than 10%. This usually requires only a 
relatively small number of points at low frequencies. 
However, if a(fc max ) <C 1, a large number of points is 
needed at high frequencies in order for the Monte Carlo 
integrator to resolve these small scale structures. For 
comparison, the instantaneous synchrotron power, Equa- 
tion (|22p. is also calculated at each point. 

Figurcs[2]-[4]show the spectra produced by an isotropic 
homogeneous particle distribution in turbulent isotropic 
fields with zero mean field component. A common fea- 
ture of each of these spectra is a hardening at low fre- 
quencies. This arises because the particle begins to be 
deflected by the turbulence through an angle compara- 
ble to that of the beaming cone of the radiation, whilst 
traversing a photon formation length, which grows to- 
wards low frequency. If the particle motion can be de- 
scribed as diffusion in the (small) angle between its ve- 
locity vector and a suitably chosen coordinate axis, this 
is kn own as the Landau-Pomeranchu k- Migdal (LPM) ef- 
fect (jLandau fe Pomeranchuklll953allbl ; iMigdali Il956l ) - 
a well-studied phenomenon in the context of the sup- 



pression of bremsstrahlung and pair-produ ction in crys - 
tals and other media (for a review see iKleml [1999I) . 
though not usually considered in the context of syn- 
chrotron radiation (alth ough see iToptvgin fc Fleishman! 
fl987t lFleishn^ar]l2006bf ). The effect can be understood 
as follows: For a trajectory with constant curvature k, 
and constant acceleration, /3 = 0, the particl e dis placc- 
mcnt is A w (3t — c 2 k 2 t 3 /24 (see equation IA8[) . For 
low frequencies the resulting formation length is dom- 
inated by the r 3 term. It is this scaling that gives 
synchrotron radiation its cj 1 / 3 asymptote at low fre- 
quencies. However, in turbulent fields, both the accel- 
eration and curvature vary. Thus, at low frequencies, 
when the formation lengths are long, a particle can un- 
dergo multiple scattering within a formation time. In 
general, for small angle scattering (a <C 7), the dis- 
placement is A w fir - \ ^dt9 2 (t) + ± dt9(t)) 2 , 
and the spectrum shou ld be averaged over a large en- 
semble of particles (e.g. | Landau fc Pomeranchuklll953at 
lAkhiezer fc Shurgalll987t) . For pitch-angle diffusion, i.e. 
(9) = and (6* 2 ) oc t. the average displacement is pro- 
portional to t 2 at low frequencies, resulting in an a; 1 / 2 
spectrum. However, the transition to this regime requires 
many scatterings and may not be realized within a for- 
mation length in a specific realization of a turbulent field. 

This is illustrated by the examples described in the 
following subsections. The radiation spectra produced in 
these examples can be placed into three broad categories, 
corresponding to the two extreme cases where a(fc m i n ) 3> 
1 or a(fc m ; n ) <C 1 and an intermediate range in which 
a(fcmin) is of order unity. 

3.4.1. a(fc min ) > 1 

In the case of large strength parameters, since the par- 
ticle in general sweeps through an angle larger than its 
beaming angle, the spectrum should resemble that of the 
instantaneous synchrotron spectrum close to the criti- 
cal frequency. Fig. [5] shows the resulting spectrum for 
two different field configurations. As expected, the spec- 
trum matches very closely that of the instantaneous syn- 
chrotron approximation in the vicinity of the roll-over 
frequency. Below this value, the numerically determined 
spectrum diverges slowly from the instantaneous syn- 
chrotron line, becoming gradually harder at lower fre- 
quencies. Note that in the large strength parameter 
regime, the transition to the diffusive LPM regime de- 
scribed above, i.e. the w 1 / 2 scaling, should occur when 
the formation length exceeds the longest wavelength in 
the system, which occurs only at very low frequencies 
uj ~ Lu c /a(k m i n ) 3 . The divergenc e fro m the synchrotron 
spectrum follows from Equation (|A8|) since now both k 

and (3 are non-zero, and the coefficient of the r 3 term will 
have an additional time-dependence. For the range of 
frequencies considered, the spectrum does not approach 
a low-frequency power-law asymptote, but continues to 
harden gradually as it approaches w g , where the syn- 
chrotron approximation fails and the beaming cone is 
large. The smaller the value of a(fc m i n ) the more rapidly 
the spectrum diverges from that of the instantaneous 
synchrotron approximation. The spectra are not sen- 
sitive to the value of a(fc max ) in the a(/c m i n ) 3> 1 regime, 
provided a(fc max ) <^ 1. 
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For frequencies above the roll-over frequency, as can 
be seen in Fig. [2l the spectrum is in excellent agreement 
with the instantaneous synchrotron approximation. On 
physical grounds it is expected that a power-law tail must 
occur at higher frequencies due to the high frequency 
jittering resulting from modes with a(k) < 1 (see e.g. 
IFleishmanl 12006a) . However, for the turbulent spectra 
considered, the power associated with such fluctuations 
is extremely small, and the numerical accuracy required 
to resolve such a feature in the large a(fc m i n ) regime is be- 
yond the capabilities of current computational resources. 

3.4.2. a(fc min ) < 1 

For fields composed exclusively of small strength pa- 
rameter fluctuations, the particle deflections are small 
and, at sufficiently high frequencies, i t is possible to use 
standard perturbation techniques (lLandau fe Lifshit3 
119711: iMedvedevi [20001: IFleishmanl l2006bO . Numerically, 
this regime is far more challenging since the time steps 
in the integrator must resolve deflections in the parti- 
cle's trajectory on the order a(fc max )/7. Figs. [3] and [4] 
show the spectra produced in fields with a(fc m i n ) = 0.5 
and a(fc m i n ) = 0.1 respectively. Both spectra exhibit a 
break at the critical frequency u> sa "f 2 k m - m c. Above this 
frequency, the spectrum has a power-law slope matching 
that of the turbulence spectrum. This can be understood 
as the up-scattering of the virtual photons of the field by 
the mono-energetic particles. In principle, the power law 
should extend up to uj rj 7 2 fc max c, however, the numeri- 
cal accuracy of the integrator chosen for this example is 
insufficient to display the entire range. At even higher 
frequencies, uj 3> 7 2 fc max c the fields are constant over the 
formation length of the particle, and the instantaneous 
synchrotron approximation applies. As we discuss in sec- 
tion 31 if the formation length of a particle is not well re- 
solved, it is exactly in this regime that the instantaneous 
synchrotron approximation must be used. Below the crit- 
ical frequency, the photon formation time remains short 
compared to the time taken to deflect through an angle 
greater than 7 . The displacement is approximately 
A ss (3t and, as in the case of relativistic bremsstrahlung, 
the spectrum is approximately flat dW/dw cx u>°. Ulti- 
mately, at frequencies u < a(k m - nl )uj c « a(fc m in) 2 7 2 cfc m in, 
the formation time exceeds the time needed to diffuse 
out of the beaming cone and the spec trum is determined 
by the LPM effect (|Fleishmanll2006bD . 

3.4.3. a(fc min )~l 

The intermediate range where the value of a(fc m ; n ) 
is somewhat larger than unity, is interesting because it 
eme rges from PIC simulations of Weibel mediated shocks 
(e.g iSironi fc~S pitkovskv 2009b). An example of the 
spectrum produced in such a field is shown in Fig. [5l For 
this example, the strength parameter a(fc m ; n ) is of order 
unity, and the transition to the LPM regime occurs at rel- 
atively high frequencies, close to the roll-over frequency. 
However, above the roll-over frequency, unlike in the 
a(^min) 3> 1 regime, the small strength parameter modes 
can be resolved, and similar to the a(fc m j n ) <C 1 spec- 
tra, a high frequency power-law emerges. The shape of 
the power-law matches that of the turbulence spectrum. 
This presents a possible observational signature of short 
wavelength turbulence at relativistic shocks. The pres- 
ence of such short wavelength turbulence is supported 



by current PIC simulations in which Fermi acceleration 
is found to occur. 

3.4.4. Non-zero mean field 

In general, the radiation spectrum can be affected if the 
mean field is non-zero or if turbulence is generated on dif- 
ferent scales such that the large scale fluctuations act as 
a local mean field. The latter situation could in principle 
be realized in the presence of large scale MHD turbulence 
produced from interaction be tween the shock front and 
density inhomogeneities (e.g. ISironi fe Goodman! [2001 
and short wavelength turbulence produced via kinetic 
effects in the shock transition region. In the presence 
of two populations of scatterers, if they are generated 
on very different length-scales, it is possible for the syn- 
chrotron radiation of shock-accelerated particles to ex- 
tend into the gamma-ray range, whereas for a single pop- 
ulation of scatterers radiation losses restr ict it to rela- 
tively low frequency (jKirk fc Re villd 12010ft . 

Here we consider the radiation produced in a region 
with a mean field having a superimposed turbulence 
spectrum. If the energy in the turbulent fluctuations is 
negligible with respect to the total field, 77 <C 1, where 77 
is defined in the low frequency spectrum will match 
that of the instantaneous synchrotron spectrum, since 
scattering will be ineffective, and to zeroth order, the 
particles simply gyrate about the mean field. At higher 
frequencies, provided modes with a < 1 exist, a power- 
law tail can emerge. Again, depending on the power 
associated with these modes, the numerical scheme can 
capture this feature, provided it is not too deep in the 
exponential cut-off region. To illustrate this, we show in 
Figure [6] the spectra produced in turbulent fields, with 
modest maximum strength parameters, and different val- 
ues of 77. For small 77, i.e. weak turbulence, the spectrum 
reproduces that of the instantaneous synchrotron spec- 
trum, although a high frequency tail is also produced, 
due to the fluctuations on modes with a < 1. As 77 in- 
creases, more power goes into the high frequency emis- 
sion, and a reduction in the power at low frequencies is 
observed, although for the frequencies investigated, the 
spectrum maintains a u; 1 / 3 scaling. 

As the ratio of the energy density in the turbulent field 
to that of the total field is increased further, we return 
to the previously investigated regi mes. For example, in 
a Weib el mediated shock 1 — 77 <C 1 (jSironi fc Spitkovskvl 
l2009a| ). However, to investigate clearly identifiable signa- 
tures, it is necessary to move beyond the prescribed, ho- 
mogeneous turbulent fields considered here to more self- 
consistent realizations, resulting from the simulations. 

4. TRAJECTORIES AND FIELDS GIVEN AS A 
TIME SERIES 

A PIC simulation is capable of producing a large num- 
ber of time series listing the position and velocity of the 
simulation particles and the values of the electromag- 
netic fields at each time-step. Using these, it is possible 
to produce spectra and light curves. The synthetic spec- 
tra presented in section [3] are based on isotropic mono- 
energetic particle distributions as described in section 
13.11 In general, however, the particle distribution is not 
only energy dependent, but can be highly anisotropic. 
This can in principal be studied by numerically solving 
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Fig. 3. — Radiation spectra for an isotropic distribution of par- 
ticles in a fully turbulent field with all strength parameters a < 1. 
The blue line is the integrated instantaneous synchrotron spectra. 
The high frequency asymptote is close to the shape of the turbu- 
lent spectrum oc oj -11 / 3 . The low frequency spectrum does not 
converge to a power law for the range of frequencies considered. 




Fig. 4. — Radiation spectra for an isotropic distribution of parti- 
cles in a fully turbulent field with strength parameters a < 1 using 
a smaller value for a(fc m ; n ) and a larger dynamic range than in 
Fig. \3\ The blue line is the integrated instantaneous synchrotron 
spectra. The high frequency asymptote is close to the shape of 
the turbulent spectrum oc o>~ 3 . The low frequency spectrum has 
a spectral slope ~ 0.7. 



Equation ([T]) or Equations (|2|) and Q for each trajec- 
tory, and then summing over trajectories, which is equiv- 
alent to integrating over the particle distribution func- 
tion. However, the radiation from an individual trajec- 
tory is beamed into an opening angle ~ If this 
is smaller than the scales on which the particle distri- 
bution is anisotr opic, the order of these operat ions can 
be reversed (see iGinzburg fc Svrovatskiil 119651 section 
3.2). The average over the particle distribution is then 
replaced by an integration over angles of the radiation 
emitted by a single trajectory (which can be performed 
analytically), and the radiation observed in a given vir- 
tual detector is given by summing over all those trajecto- 
ries whose velocity vector lies within the acceptance cone 
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Fig. 5. — Radiation spectra for an isotropic distribution of parti- 
cles in a fully turbulent field with strength parameter a(fc m in) > 1 
and a(fc max ) <§C 1. The blue line is the integrated instantaneous 
synchrotron spectra. The index of the high frequency power-law 
component is close to that of the turbulent spectrum oc a; -8 / 3 . Ev- 
idence of a cut-off is observed close to where the formation length 
L c ~ l/k ma x, where the line must match up with the instanta- 
neous synchrotron approximation. This cannot be resolved numer- 
ically. 
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Fig. 6. — Spectra emitted in the presence of a finite mean field 
r\ = 0.1 (red curves) and 7] = 0.9 (blue curves), where r\ is the ra- 
tio of the ene rgy density in the turbulent field to the total energy 
density — see l|35J l. The instantaneous synchrotron spectra are plot- 
ted using dashed lines. The total energy density in the magnetic 
field is fixed, so that the magnitude of the average field differs in 
the two cases. As in Fig. [5] there is evidence of a cut-off at high 
frequencies. 



of that detector. Formally, 
dL 



dudn 



= J d 3 xdpp 2 d 2 flf{x, P fl,t)P(n,u},t) 
w J d 3 xdpp 2 f(x,pn,t) J d 2 ClP(Cl,uj 7 t) 

= J d 3 xdpp 2 f(x,pn,t)P(u,t) (38) 

and the integrations over x and p reduce in the PIC case 
to summations over all trajectories that illuminate the 
specified detector. 
For the high-energy emission of particles accelerated at 
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a rclativistic shock front, the restriction imposed by this 
procedure is not important, because the anisotropy of the 
particle distribution is expected to be on a scale larger 
than the beaming angle. Thus, the angular dependence 
of the emitt ed radiation foun d bv [ Sironi fc Spitkovskyl 
(|2009bft and lFrederiksen et al.1 ((20101) should just reflect 
the angular dependence of the distribution function at 
the relevant particle energy, and would be preserved in 
this approach. 

As pointed out by iHededal (|2005l ) the computation of 
synthetic spectra from trajectories taken from PIC simu- 
lations inevitably involves interpolation. Specifically, the 
algorithm presented in (f9"|- (fTTj) transforms the integra- 
tion variable from time to phase. In order to split the 
contributions to the integral into an alternating series 
([2"5")h the discrete trajectory must be interpolated. 

Interpolation is not a sensitive procedure provided 
many points are contained within a photon formation 
time, a constraint that will be made more precise be- 
low. An accurate evaluation of the instantaneous power 
at any time step can, for example, be obtained simply by 
linearly interpolating the functions <?, 7, (3, g, 8j3 and 5 A, 
which arc known at all neighboring grid points. When 
the photon formation length drops to only a few time 
steps, this procedure fails. However, the validity of the 
PIC simulation requires that the electromagnetic fields 
vary slowly between time steps, which is precisely the 
condition for applicability of the generalized synchrotron 
formula (|22|) . Therefore, in a valid simulation, the instan- 
taneous power can safely be evaluated using this method, 
if the formation time is not long compared to the time 
step. It follows that, for a given frequency, the method 
of evaluating the instantaneous power at each of the dis- 
crete set of particle positions x(t n ), depends on the value 
of the photon formation time at that point. 

At high frequencies, the formation time is short, and 
can be much shorter than the typical time-step used in 
PIC simulations, which is a fraction of a plasma cycle. It 
is straightforward to find for each time-step (labeled by 
n) the values 8A^ of the deviation of the displacement 
at the neighboring points nil. For a given frequency, 
the photon formation lengths in the forward and back- 
ward directions follow. Alternatively, two critical fre- 
quencies can be found such that at these frequencies 
the neighboring points lie precisely one formation length 
away from x n . From Equations (|4]) and (fT4|) the critical 
frequencies are 
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2 7 2 |<5A, 



(39) 



where St^ = \t n ±i — t n \ is the time-step between neigh- 
boring data points. If u> is close to or greater than oj^, 
then the coherence length is poorly resolved and the syn- 
chrotron approximation must be used to compute the in- 
stantaneous power. If, on the other hand, ui <C w^, then 
the coherence length is well resolved, and a numerical 
integration is accurate. 

The accuracy of the solution depends quite strongly on 
the ability to resolve the peaks and troughs of the sine 
function in Equation (|25p . PIC simulations usually work 
with a fixed time-step, in which case the resolution in 
successive terms in ([25)) decreases. This case be seen by 
considering the time evolution of the phase. Making a 



Taylor expansion about the initial position gives 



g « u 



J_ + -1 c 2 /3 3 k 2 t 3 
2 7 2 24 H 



(40) 



For t > 3/7CK, the r 3 term dominates and one can solve 
for g = mr to give 
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(41) 



It is readily seen that for larger n the time interval be- 
tween successive integer multiples of the phase g = nir 
decreases: 



'ra+l 



T 



2tt lo c 



1/3 



(42) 



For this reason, it is essential to interpolate the functions 
Q and g, rather than the combination Qsing, and we 
have found that linear interpolation is adequate. Then, 
trapezoidal integration in the phase g is used with a max- 
imum step-size of Ag = 27r/25fl to evaluate the terms 
n= -10... 10 in (©.' 

As an illustrative example, we again consider the case 
of uniform circular motion. In Fig [7] we plot the energy 
radiated per frequency interval over one gyration. In ad- 
dition to the analytic solution, the result of integrating 
the instantaneous power calculated using various time- 
steps is shown. In this special example, both the instan- 
taneous power and the frequencies uj^ are independent of 
time, so that use of the synchrotron approximation au- 
tomatically yields the exact analytic result. The numer- 
ically determined power reproduces this result to within 
1% for frequencies 
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(43) 



At higher frequencies, the instantaneous power itself 
may still be evaluated accurately, since the interpola- 
tion scheme guarantees 25 points per photon formation 
length. However, since this quantity is evaluated only 
at each time-step, the subsequent integration required to 
evaluate the radiated energy does not reach the required 
resolution. 

This suggests the following procedure when the algo- 
rithm is employed in an arbitrary field configuration: The 
frequency at which the emission is to be evaluated, is 
compared at each time-step to the frequencies oj^. If 
u satisfies the inequality (|43|) . numerical integration is 
used. Otherwise, the instantaneous synchrotron expres- 
sion (|2"2")) is used. 

An important property of this algorithm is that it 
avoids explicitly interpolating the particle's position and 
velocity. Such a procedure introduces discontinuities into 
the particle acceleration as a function of time, leading to 
high-frequency artifacts similar to those that arise when 
the accelera tion of a hyperbolic trajectory is abruptly 
terminated (jReville fc Kirkll2010[) . 

As mentioned above, the instantaneous power can only 
be evaluated accurately by integrating over at least the 

1 This is approximately the resolution required to calculate 
f a s'm(x)dx to better than 99% accuracy using trapezoidal inte- 
gration. 
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Fig. 7. — Synchrotron spectra found using discrete time-series 
data for a 7 = 10 3 particle using linear interpolation on Q and 
g between data points, for a range of time-step sizes (in units of 
ujg ). In each case, the frequency at which the numerical result 
begi ns to deviate from the exact answer is in close agreement with 



first few formation lengths. The number of terms needed 
in (|25|) can be considerably reduced with the aid of the 
Euler-van Wijngaarden transform. However, when cal- 
culating the total energy spectrum radiated by an indi- 
vidual trajectory, it is also essential to resolve the instan- 
taneous power as a function of time. Given a finite time 
series of positions and velocities, the radiation formulas 
apply only if the trajectory is extrapolated ballistically 
outside of the finite length time series, although the be- 
havior in these regions does not affect the results when 
frequencies cj 3> (-f 2 )/T are considered. Here (j 2 ) is the 
average Lorentz factor squared along the trajectory, and 
T the total time. Of course, if a time-dependent light- 
curve is to be generated, the restriction is much more 
severe, since then T refers to the time-interval between 
successive evaluations. 

To illustrate this, we consider the spectrum produced 
by a rclativistic particle that undergoes an instantaneous 
scattering at t = through an angle a. The spectrum 
in this case is well known to be flat, u , at frequen- 
cies small compared to the inverse du ration of the ac- 
celera tion (for a detailed discussion see iSchwinger et al.l 
119981 chapter 37). The angular integrated spectrum in 
this frequency range can b e determined analytically (e.g. 
lAkhiezer fc Shul'galH98l 
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Although this example appears straightforward, it is, 
in fact, quite demanding numerically, the reason for this 
being that the instantaneous power is itself an oscilla- 
tory function. The formation length at any given time t 
can be easily calculated, but the exact expression is cum- 
bersome. Far from the scattering center, t c = Air-f 2 /u. 
As the scatterer is approached, the formation length de- 
creases, reaching a minimum at t = 0, of 
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Fig. 8. — Instantaneous power at frequency u> = 1 (in arbitrary 
dimensionless units) as a function of time produced by a particle 
with 7 = 10 3 that undergoes an instantaneous scattering at t = 
through an angle a = 2 X 10 — 4 (£ = 0. 1). Time is measured in 
units of the coherence time given by Eq. H45II . 



In this example there is no intrinsic time scale, so that we 
are free to choose arbitrary time, distance and frequency 
units. Defining a reference time unit t$, we construct 
dimensionless units t = t/to, x = x/cto and Co = uto. 

Figures [8] and [9] show the instantaneous power as a 
function of time for £ = 0.1 and £ = 10, respectively. The 
power oscillates with a slowly increasing period approx- 
imately equal to the formation time and damping with 
distance from the scatterer. In addition, there is an un- 
resolved discontinuity at t = 0, which arises because the 
particle velocity is also discontinuous at this point. How- 
ever, after integration over t, this feature has no influence 
on the energy radiated. Clearly, the linear growth phase 
for Q{g,t) will increase with distance from the scatter- 
ing event, and the number of terms in the Euler-van 
Wijngaarden transform should be chosen such that the 
scattering is included. However, since the dominant con- 
tribution to the total energy radiated comes from the 
first few periods, the integral of the instantaneous power 
converges rapidly and the error incurred from taking only 
the first few formation lengths when calculating P(io,t) 
is small. 

To demonstrate the effects of having a finite trajectory, 
we integrate the instantaneous power over a time interval 



dE 
duj 



P(io,t)dt 



—T 



(45) 



with T = 7T7 2 /(1 + 4£ 2 ), corresponding to one forma- 
tion length for an emitted wave with frequency u = 4. 
From figures [8] and [9l it is clear that the solution will 
converge only for frequencies much larger than this. The 
instantaneous power is integrated using a finite time step 
trapezoidal integration for two different scattering angles 
£ = 0.1 and £ = 10, with 7 = 10 3 , as above. The results 
are shown in figure QUI The result is in good agreement 
with the analytic solution above approximately u = 20. 
This suggests that for a given trajectory, on a time inter- 
val [— T, T], the minimum frequency that can be inves- 
tigated must have at least 10 formation lengths in this 
time interval. 
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Fig. 9. — Instantaneous power at frequency cj = 1 as a function 
of time produced by a particle with 7 = 10 3 that undergoes an 
instantaneous scattering at t = through an angle a = 2 X 10 — 2 
(5 = 10). Time is measured in units of the coherence time given 
by Eq. @5) . 




Fig. 10. — Plot demonstrating the low frequency limitation due 
to finite end poin ts of the trajectory. The dashed line is the analytic 
result from 144 p . For the numerical evaluation, the instantaneous 
power was evaluated on the time interval —T<t<T where 
T = 7T7 2 /(1 + 4£ 2 ) with 7 = 10 3 and £ = 0.1 and £ = 10. 



5. DISCUSSION 

In this paper we describe an algorithm for calculating 
the radiation emitted by a relativistic charged particle 
moving in turbulent electromagnetic fields, and use it to 
investigate the spectra that arise in a prescribed, stochas- 
tic realization of a static, turbulent magnetic field. Wc 
also describe how to adapt the approach when the trajec- 
tory is given as finite time series for the position, velocity 
and acceleration. The algorithm is based on formulating 
an "instantaneous power" at each point on the trajec- 
tory, and makes use of the concept of the photon forma- 
tion length to evaluate this quantity. It is suitable for 
use in post-processing the output from a particle-in-cell 
simulation. 

Two problems arise with a numerical evaluation of the 
radiation. At high frequencies, the finite time-resolution 
of the trajectory is a limitation. For relativistic particles, 
this problem can be alleviated by suitable grouping of 
the terms associated with the slowly and rapidly varying 



components of the instantaneous power, which improves 
the stability and accuracy. But even with the appro- 
priate grouping of the terms, if the time-resolution of 
the trajectory cannot be improved indefinitely, a purely 
numerical evaluation still fails at sufficiently high fre- 
quency. Fortunately, it is precisely in this range that the 
instantaneous power can safely be evaluated using the 
synchrotron approximation. Subsequent integration of 
this quantity does not present a difficulty. However, for 
non-relativistic particles, or for frequencies comparable 
to the instantaneous angular frequency of the emitting 
particle, additional terms enter into the expression for 
the instantaneous power (|15|) 

At low frequencies, a limitation is imposed when the 
trajectory is known only within a finite time interval. 
This intrinsic restriction cannot be removed: the mini- 
mum frequency at which a light curve can be computed 
is roughly I0j 2 /T, where T is the length of the available 
time series. In current PIC simulations, however, the 
neglect of the collective response of the plasma to the 
propagating waves, which leads to effects such as Razin- 
Tsytovich suppression and transition radiation, is also 
likely to be important. These should intervene at fre- 
quencies below roughly ~ juj p , where ui p = y/ 'Airne 2 / 'm 
is the plasma frequency and n the number density, but 
to date it does not appear feasible to account for such 
effects self-consistently. For particles of Lorentz factor 
100, these estimates imply that the intrinsic restriction 
is more important than the neglect of collective effects 
only when time-series shorter than about 10 3 plasma cy- 
cles are used to compute the emitted radiation. 

Our results obtained for 3D magnctostatic turbulence 
confirm that one observational signature of short length- 
scale turbulence, is the presence of a high-frequency 
power-law tail in the mono-energetic emission spectrum 
iFleishmanl ()2006af ) . For a power law of electrons dn/d7 cx 
7~ p , one expects a synchrotron power-law spectrum 
F u oc uj~ s , where s = (p — l)/2, for frequencies be- 
low the roll-over frequency of the maximum energy elec- 
trons. Observations of GRBs place this index in the 
range 2 < p < 2.8, corresponding to a spectral index 
of 0.5 < s < 0.9. Thus, unless the turbulence index is 
extremely hard, a < 1, the photon spectrum will not 
harden at high frequencies, and this observational signa- 
ture may be difficult to distinguish from a cut-off. 

The synchrotron spectrum of particles radiating in 
a uniform field is nowhere harder than an w 1 / 3 power 
law, which, since harder spectra have been observed in 
gamma-ray bursts, has l ed to the discussion of a syn- 
chrotron 'line of death' (|Preece et al.l [199 8) . Our re- 
sults confirm that this is generally true for isotropic 
particle distributions in large scale, static, 3D turbu- 
lence. However, in agree ment with other treatments 
(|Fleishman fc UrtievlfeOlOj) we find that, in the presence 
of large amplitude turbulence on short length-scales, the 
low-frequency asymptote can diverge from this value. We 
find low frequency spectra that are typically harder than 
w 1 / 3 . For fields with a 1 the spectrum exhibits a grad- 
ual hardening with the slope increasing approximately 
0.05 per decade in frequency. For fields with strength 
parameters a > 1, low frequency power-law asymptotes 
are produced, F u oc uj q , with 1/3 < q < 1/2. For fields 
with a < 1 slightly larger values of q appear to be possi- 
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ble, although shocks w ith such small strengt h parameters 
arc poor accelerators (|Kirk fc Revillell2~010[ ). Spectra as 
hard as oj 1 are known to be produced by weak (a <C 1) 
turbulen ce that can be fac t orized into 2D a nd ID com- 
ponents ([Fleishman! I2006tj ; iMedvedevl 120061) . However, 
they do not arise in our results, which are based on a 



fully 3D turbulence model. 



We thank A. M. Taylor and S. O'Sullivan for helpful 
discussions. B.R. gratefully acknowledges support from 
the Alexander von Humboldt foundation. 
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APPENDIX 



SYNCHROTRON EMISSION 

We start by making a Taylor expansion of the particle position, the deviations from ballistic motion and the phase. 
This can be achieved using a purely geometric description of the trajectory. 
First, define the tangent, normal and binomial unit- vectors: 

AT /At AT /At 



N ®' \AT/At\ cP(t) K (t) 
B(t) — T(t) A N(t) (Al) 



The quantity K(t) is called the curvature of the trajectory. The Frenet-Serret formulae give the evolution of these 
vectors along the trajectory: 



/dT/dt\ 


/ c/3k \ 


n 


AN /At -- 


—c/3k c/3t 


N 


K AB/At ) 


V -c/3f / 


\b) 



where r is called the torsion of the trajectory. 



Therefore, 



Then, using 



SR(t, t) 



r r 2 r 

—c 2 P 2 k + — { 3c 2 f3/3n + c 2 (3 2 k 



N+ 1T C ' 
6 



r/3 + ^-cW) 



c z f3 6 n?B + 0(r d ) 



T + 



tc/3 2 k + — [3c/3(3k + c/3 2 k 



« 3 KfB + 0(r 4 ) 
N 



2 3 

(3 ■ 5R/c= T -pl3 + T - (pp - c 2 /3V) + 0(r 4 ) 



(&Rf /<? = ?- 



(/3 2 + c 2 /3 4 K 2 ) + 0(r 5 ) 
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(A3) 

(A4) 

(A5) 
(A6) 



one finds 



5A(t, r) = \t 2 /3 2 + 2t(3 ■ SR/c + (SR) 2 j c 2 ] 1 ' 2 - 
= ^ + M(4/3- c W)+0(t 4 ) 



Substituting into the definition of the phase-lag 

g(t,r)=UT 



T 

24 



^K 2 



g = LO 



+ 0(r 4 ) 
+ 0(r 3 ) 



(A7) 
(A8) 

(A9) 
(A10) 



U$ - c 2 /3 3 

l-P-rp- T - ^/3-c 2 (3 3 k 2 ' 

At this point, two additional assumptions are introduced: 

1. the electromagnetic fields are constant over a photon formation length, i.e., ft = 

2. linear acceleration emission is negligible, i.e., /3 = 

The first is an implicit condition for the validity of a PIC simulation, when the photon formation time is comparable 
or shorter than the time step. The second is fulfilled under normal conditions (\E\ < \B\). 
Then, writing 



CKJT 



-7 CK 



(All) 
(A12) 



one finds 



g(t,r) ~ujt 

3w 
2u> c 



1 ~ /3+ 2i° P KV 



and 



(A13) 
(A14) 

(A15) 



At the frequencies of interest (cj ~ w c ), the dominant contribution to the integrals, which arises for g ~ 1, occurs for 
2: ~ 1. In this case, the the contribution of the first two non- vanishing terms in the Taylor expansions of both g and 
g are comparable, Thus, in expanding the integrands in (|18p it is necessary to include both these terms, whereas the 
lowest order non- vanishing contributions to <5A and f3 ■ 8(3 arc sufficient. 
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The instantaneous power (|18[) is then 



Writing rj = lu/uj c , 



Then, using 



to find 



and noting that 



e 2 uj 4 



oo x 



TIC 7- J 



dx 



„2 , x" 



3oj ( x A 



dx < ~ 2 



cos 



2 V 3 



dx x sin 



377 / x A 



dx x 2 cos 



da; x cos 



377 / x s 
— [ x H 

3l] ( x 3 



1 

2 



— ^3(7?) -3X 2 



d , . 2 



one arrives at the standard expression for angle-integrated synchrotron radiation 

P(V,t)= 2n V dxK 5/3 (x) 



(A16) 

(A17) 
(A18) 

(A19) 
(A20) 

(A21) 
(A22) 



